FINITE ELEMENT ANALYSIS OF THE FORCED 
OSCILLATIONS OF SHIP HULL FORMS. 



Davfd Arthur Smith 



X library 

iRADUATE SCHOOL 
CALIFORNIA 93940 




li^onterey, California 





FINITE ELEMENT ANALYSIS OP THE FORCED 
OSCILLATION OP SHIP HULL FORMS 

by 

David Arthur Smith, Jr. 

June 197^ 



Thesis Advisor: 



R.E. Newton 



Approved for public release; distribution unlimited. 



security classification of this page Dmtm Entt*rt*d) 



REPORT DOCUMENTATION PAGE 


READ INSTRUCTIONS 
BEFORE COMPLETING FORM 


1. REPORT NUMBER 


2. GOVT ACCESSION NO. 


3. RECIPIENT'S CATALOG NUMBER 


4. Tl TL E (’and Su6rif/«; 

Finite Element Analysis of the Forced 
Oscillations of Ship Hull Forms 


5. TYPE OF REPORT 4 PERIOD COVERED 

Engineer's Thesis; 

June 197^ 


6. PERFORMING ORG. REPORT NUMBER 


7. AUTHORf*; 

David Arthur Smith, Jr. 


B. CONTRACT OR GRANT NUMBERf*) 


9. PERFORMING ORGANIZATION NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, California 939^0 


10. PROGRAM ELEMENT. PROJECT, TASK 
AREA & WORK UNIT NUMBERS 


II. CONTROLLING OFFICE NAME AND ADDRESS 

Naval Postgraduate School 
Monterey, California 939^0 


12. REPORT DATE 

June 197^ 


IJ. NUMBER OF PAGES 
102 


14. monitoring agency name & ADDRESS^// different from Controlling Office) 

Naval Postgraduate School 
Monterey, California 939^0 


IS. SECURITY CLASS, (ol Ihlm fSjwrO 

Unclassified 


15a. DECLASSI FI cation/ down GRADING 
SCHEDULE 



16. DISTRIBUTION STATEMENT (ot thia Report) 



Approved for public release; distribution unlimited. 



17. DISTRIBUTION STATEMENT (o( the mbetract entered in Block 20, H different from Report) 



18. SUPPLEMENTARY NOTES 



19. KEY WORDS (Con^/nua on reverae aide If necaaemry end identify by block number) 

Added mass 
Finite element 
Nonlinear 



20. ABSTRACT (Continue on reverae aide If neceaaery and Identify by block number) 

A study of the first and second-order, dynamic, ideal fluid 
effects on the forced, harmonic oscillations of a rigid cylinder 
in a free surface is presented. The solutions are obtained v;ith 
a FORTRAN IV computer program based on the finite element method 
using Isoparametric elements. 

First-order added mass and damping in heave for a semi-circular 
and bulb hull form including finite depth effects are presented. 
Second-order solutions for the semi-circular hull form in heave 



DD 1 473 

(Page 1) 



EDITION OF 1 NOV 65 IS OBSOLETE 
S/N 0102-014- 6601 | 



1 



SECURITY CLASSIFICATION OF THIS PAGE (Wxtn Dmtm Bntm 



ClCUKITY classification of This PAOErWTien Dmtm Enffd) 



(20. ABSTRACT continued) 

also including finite depth effects are obtained. The solutions 
are verified by comparing with existing theory and experiment. 



DD Form 1473 
, 1 Jan 73 

S/N 0102-014-6601 



2 



SECURITY CLASSIFICATION OF THIS PAGEfWTion Dmtm Entmrmd) 



Finite Element Analysis of the Forced 
Oscillations of Ship Hull Forms 

by 

David Arthur^^Smith , Jr. 
Lieutenant, United States Navy 
B.S., University of Missouri, 1966 



Submitted in partial fulfillment of the 
requirements for the degrees of 

MASTER OF SCIENCE IN MECHANICAL ENGINEERING 

and 

MECHANICAL ENGINEER 



from the 

NAVAL POSTGRADUATE SCHOOL 
June 197^ 









t 



ABSTRACT 



I 



A study of the first and second-order, dynamic. Ideal 
fluid effects on the forced, harmonic oscillations of a 
rigid cylinder In a free surface Is presented. The solutions 
are obtained with a FORTRAN IV computer program based on the 
finite element method using Isoparametric elements. 

First-order added mass and damping In heave for a semi- 
circular and bulb hull form Including finite depth effects 
are presented. Second-order solutions for the seml-clrcular 
hull form In heave also Including finite depth effects are 
obtained. The solutions are verified by comparing with 
existing theory and experiment. 
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I. INTRODUCTION 



Naval architects and marine engineers have long been 
greatly interested in the study of the dynamic fluid response 
to the forced oscillations of a two-dimensional floating body. 
The applications of the studies include studying the rigid- 
body response of the ship due to an incoming wave train, and 
determining the natural frequency and vibration mode shapes 
of a ship's structure in connection with propulsion systems 
design [3^1* 

Marine engineers, in modern times, are concerned with, 
among other things, the natural frequencies and stability 
characteristics of floating platforms near a coastline. 

Also, they are concerned with the related problem of finding 
the dynamic loads on submerged, moored, or fixed bodies [8]. 

The modern techniques for solving these types of problems 
involve seeking the solution of a steady-state, periodic, 
potential flow problem with a free surface in a fluid of 
either infinite or finite depth. This problem is made 
tractable by linearizing the free surface boundary condition. 

A. HISTORICAL BACKGROUND 

The historical development of potential flow problems 
involving a free surface with wave motion, and the related 
problem of the dynamic response of fluids due to the motion 
of rigid bodies on and under a free surface is lengthy; but 
of interest. Several authors have reviewed this history in 
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some detail [3,17,2^1,31]. The history is repeated in the 
following in the interest of providing a foundation upon 
which this author endeavors to add additional Information in 
an attempt to provide a more complete history. 

1 . First Order Solutions 

a. General Potential Theory 

The foundation of the modern day theory of fluid 
motions begins, of course, with the classical works of Newton, 
Lagrange, Euler, Stokes, and others. 

However, in 1879, Lamb [15] provided one of the 
first complete studies of the theory of wave motions; which 
he updated five times until 1932. It appears that the next 
major contribution was made by Havelock [10] in 1929; wherein, 
he studied the solution of the wave forms created by a wave 
generator . 

The modern beginnings of the study of the inter- 
action of a fluid and a harmonically oscillating rigid body 
at the free surface start with the work of Lewis [19] in 1929. 
Lewis studied the determination of added mass with a boundary 
condition requiring zero pressure at the mean position of the 
free surface, thus eliminating the effects of frequency of 
body oscillation and of damping. Lewis considered the high 
frequency case of the more general problem. 

The first work which Included the effects of 
frequency and damping was done by Ursell [38] in 19^9. This 
work solved the problem of a right-circular, seml-lmmersed 
cylinder heaving on the free surface of a fluid infinite in 
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extent. Ursell presented added mass and damping coefficients 
as a function of frequency. His method of solution involved 
placing an "infinite" number of singularities of all orders 
at the Intersection of the free surface at rest and the 
cylinder centerline. The strengths of the sources were 
determined from the kinematic boundary condition applied on 
the Immersed surface of the cylinder. Grim [9] arrived at a 
similar solution in 1953. 

A period of about ten years passed before any 
additional work was reported on the subject. There is a 
rather obvious reason for this time lag. Ursell 's solution 
involved laborious numerical calculations which had to be 
done by hand. The advent of the digital computer made it 
feasible to pursue further solutions . 

In i 960 , Tasal [ 36 ] extended the work of Ursell 
to include roll and sway motions of Lewis forms. At the same 
time. Porter [30] compared linearized pressure calculations 
with experimental results, and extended the range of solutions 
to Include other hull forms. 

It appears that Yu and Ursell [^2] provided the 
first theoretical study of the effect of finite depth on 
two-dimensional solutions in 196 I. This was followed by 
additional work by Kim [1^] in 1969> who also Included finite 
depth. In the first work, added mass coefficients and the 
properties of the resultant linear wave were reported. The 
second work, by Kim, gave added mass and damping coefficients 
as a function of frequency and depth. Work by Paulllng and 
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Richardson [28] and Paulling and Porter [29] provided some 
experimental verification of the theory in I 962 . Vugts [^0], 
at the Netherlands Research Center in 1968, reported additional 
extensive, experimental research on seven different ship hull 
forms . 

b. Finite Element Solutions 

Zienkiewlcz [^6], in 1965> provided the first 
successful application of the finite element method to Poisson's 
equation. Although the finite element method was originally 
developed to solve problems in the theory of elastic mediums, 
Zienkiewlcz • s work "opened the door" to possible solutions to 
boundary value problems of many different types, including the 
solution of fluid flow problems; this follows since Laplace's 
equation is contained in Poisson's equation. 

In the same year (1965), Zienkiewlcz, Irons, and 
Nath [^4] first used the finite element method to find added 
nfass . Just as in the case of Lewis and others before, surface 
waves were excluded from this solution, as well as in subse- 
quent works by Rdren [32], Holand [12], Matsuura and Kawakami 
[22], and Matsumoto [21]. 

Zienkiewlcz and Newton [43] presented the first 
general theoretical treatment of the fluid-structure inter- 
action problem using the finite element theory in I 969 . Their 
work provided a means of including surface waves and deter- 
mination of energy loss due to wave propagation. The latter 
was accomplished by providing a radiation boundary condition. 
This radiation boundary condition was essential to the 
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successful modeling of Infinite regions under steady-state 
conditions. In addition, compressibility effects could be 
included in the analysis. 

Chenault [^], in 1970 , appears to be the first to 
apply the finite element analysis developed by Zienkiewicz 
and Newton to find added mass and damping coefficients as a 
function of frequency and hull form. The fluid regions he 
studied were chosen to simulate infinite depth. Some of his 
results are published in Ref. [23]. Bal [3] also Included 
surface waves in his 1972 analysis of a wide range of problems. 
His work demonstrated the flexibility of the method by studying, 
in addition to the problems of Chenault, the additions of sway 
and roll motions, finite depth, irregular fluid bottom geometry, 
axisymmetrlc geometry, and fluid stratification. Bal also 
provided criteria for properly placing the radiation boundary 
and developed a new form of the radiation condition for the 
axi-symmetric case. 

The author has extended the studies of Chenault 
and Bai in the solution of the first-order problems with an 
improved version of Chenault *s computer program. The results 
are presented later in this work along with results separately 
reported by Newton, Chenault, and Smith [23]. 

2 . Second Order Solutions 

a. General Potential Theory 

The origins of the second-order theory of gravity 
waves begin again with the work of Stokes and his classic 
solution of nonlinear traveling waves. Others have extended 
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the work of Stokes, but it appears that a modern version 
of Stokes' theory was first compiled by Stoker [35] in 1957. 
The generally accepted terminology associated with Stoker's 
work is "perturbation theory" which of course has applications 
in many fields. 

Fontannet [5], in 1961, appears to be the first 
to attempt to use the basic nonlinear perturbation theory as 
described by Stoker to solve the problem of finding the 
characteristics of waves generated by a plane wave-maker using 
Lagrangian coordinates. Ogllvie [25], in 1963, found the 
first and second-order forces on a fixed or oscillating 
circular cylinder submerged in a fluid with a free surface . 
However, his second-order forces were only those due to the 
first-order solution, and only the time-independent part was 
reported. Tuck [37] found the second-order forces on a 
subkierged cylinder in a uniform stream in 1965. Salvesen [33] 
extended Tuck's work to include wiiig-shaped bodies in his 
Ph.D. dissertation in 1966. 

Lee [17] and Parissis [27] appear to be the first 
to apply second-order perturbation theory to attempt to solve 
the fluid-structure interaction potential flow problem with 
a free surface in 1967. Lee confined his studies to the 
heaving of floating circular and U-shaped sections, while 
Parissis studied only the circular shape in heave; both 
studies were in fluids of infinite extent. These investiga- 
tions, as well as those previously mentioned, required the 
development of a fluid-structure second-order boundary 
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condition (discussed In more detail later) In addition to 
the free-surface equations compiled by Stoker. 

The work of Lee and Parissis was followed by that 
of Potash [31] In 1970. Potash extended the second-order 
studies to include the additional degrees of freedom of roll 
and sway and the attendant coupling between modes of motion. 
Recently, Garrison [7] completed a second-order solution of 
the problem of determining the dynamic loads experienced by 
a "fixed” three-dimensional body due to an incoming wave 
train. His theory Includes submerged and surface-piercing 
bodies . 

All of the second-order solutions previously 
mentioned used some form of singularity distribution in the 
flow field. The functional forms of the singular functions 
used were usually taken from the work of Wehausen and Laltone 
[^ 1 ]. 

b. Finite Element Solutions 

The author has been able to find only one 
solution to the second-order, of the class of problems being 
discussed, by the finite element method. Allouard and Coudert 
[2] presented a paper at the International Symposium on Finite 
Element Methods in Flow Problems in January of 197^. They 
presented a method of solving problems related to floating 
bodies in water waves under unsteady-linear or stationary- 
nonlinear (second-order) conditions. 
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B. OBJECTIVES OF PRESENT RESEARCH 



The experimental studies of Paulllng and Porter [29], 
Paulllng and Richardson [28], and Vugts [^0], previously 
mentioned. Indicated good agreement with the developed linear 
theory if the amplitude of motion was small. Potash [31] 
observed that small variations in amplitude caused large 
discrepancies between Vugt’s experimental results and those 
of linear theory. Potash [31] has shown that second-order 
effects become significant at higher frequencies and are more 
significant when roll and/or sway modes are present. However, 
he did not examine the effect of finite depth on second-order 
solutions . 

The present work has two primary goals. The first is to 
determine if the finite element method can be successfully 
applied to finding solutions to second-order boundary value 
problems of the type previously described. And, if success 
is achieved in the first objective, the second major goal is 
to study the significance of nonlinear effects in heave in 
water of finite depth. 
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II. FORMULATION OF THE MATHEMATICAL MODEL 



The physical problem Is to determine the force required 
to sustain steady, vertical harmonic motion of a floating 
cylindrlc ship hull form. The ship will be considered 
oscillating in an "ideal" fluid of infinite horizontal extent 
and arbitrary constant depth. 




FIGURE 1. Hull at Free Surface 

A. THE NONLINEAR BOUNDARY VALUE PROBLEM 

The domain of the boundary value problem is shown in 
figure 1. In the figure, b is the ship’s half beam, h is the 
mean fluid depth, d is the ship’s draft, and 2A is the total 
submerged cross-sectional area of the ship. 

The fluid is assumed to be Incompressible, inviscid, 
homogeneous, and without surface tension. It is well known 
that under the assumption of irrotatlonality , a velocity 
potential exists. Furthermore, the potential will be a 



single-valued function in any simply-connected region. The 
convention for the potential used herein is 



= ul + vj , 



( 1 ) 



where u is the x-component of fluid velocity, v is the 
y-component of fluid velocity, and i and are the unit 
vectors in the positive right-handed cartesian coordinate 
system directions x and y respectively. The symbol ^ is the 
standard "gradient" used in vector calculus. 

The nature of the problem dictates that at a large 
horizontal distance from the ship's hull (cross-hatched area., 
fig. 1) only an outgoing one-dimensional traveling wave will 
be present. Therefore, it is advantageous to consider the 
reduced region shown in figure 2. In figure 2, w is the 
region semi-width. The meanings of the other symbols shown 
in figure 2 will become apparent in the following discussion. 




w 






1 



Free Surface, Sj i 

^Ship Hull.Sg j 

I 



— Plane of 



Symmetry, S| ^Bottom, Sg 






FIGURE 2. Reduced Region 



Radiation 
Boundary, $4 
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Conservation of mass requires: 



= 0 , in R 



(2) 



2 

R is the fluid domain and V is the Laplacian operator. 



There are two nonlinear free surface conditions, one dynamic 
and one kinematic. The dynamic condition, from the Bernoull 
equation, is 



Equation 3 must hold on the moving free surface (p = 0) . In 
the equation, g is the acceleration of gravity, p is the 
overpressure (referred to atmospheric), P is the fluid 
density, and t is real time. The subscripts denote partial 
differentiation with respect to the indicated variable. The 
kinematic condition comes from a basic assumption from 
continuum mechanics simply stated by Stoker: "A particle 

once on the free surface remains on it." This assumption 
translates into the following expression 



which may appear, is assumed to be identically zero without 
loss of generality (See Stoker [35], p. 10). 



+ 1 [$2 + $2^ + I + yg = 0 . 




D[n(x,t) - y] ^ 
Dt 



0 



(^0 



1 



In equation 3, the possible pure function of time 
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Again, equation must be satisfied on the moving free 
surface. In equation 4, D/Dt denotes the co-movlng derivative 
(sometimes called the substantial or total derivative) and the 
function ri(x,t) Is defined by the following expression 

y = n(x,t) . (5) 

Therefore, n is defined to be the wave height above the line 

y = 0 at any point x and any time t . 

Equation can be rewritten using the definition of $ in 

2 

equation 1 and the definition of the co-moving derivative 

Vx - ty + = 0 , (6) 

on the moving free surface . 

On and "near" the bounding surface (radiation boundary) 

$(x,y,t) $*(x,y,t) , as x-»-oo . (7) 

In equation 7> denotes the velocity potential of a 
traveling wave moving away from the ship. It must be mentioned 
here that inherent in the statement of equation 7 is the 
assumption that the surface is "far enough" away from the 
ship's hull. This point will be discussed in more detail 
later . 

Turning to the interface S 2 between the fluid and the 
ship's hull, the kinematic condition is 



2 

The co-moving derivative is defined as: 
D( )/Dt E u( )^ + v( )y + ( )^ . 
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( 8 ) 



V$.n = V*n 

The vector n is the unit outward normal vector (from the 
region R) to the surface at the point of interest. The dot 
is the standard scalar product and V represents the vector 
velocity of the ship at any instant of time. 

Since the motion of the ship is assumed to be restricted 
to the vertical direction (heave), the bounding surface S 2 
becomes a plane of symmetry. The surface S^, represents the 
bottom of the fluid region and is rigid. Therefore, both 
surfaces have the same kinematic boundary condition 

^$.n = 0 . (9) 

An examination of equations 3 and 6 indicates the nonlineari- 
ties inherent in the boundary value problem. An additional 
problem, not as obvious, is the fact that the free surface 
and the hull are moving boundaries . In its present form the 
problem has not yielded solutions. Therefore, some approxi- 
mate theory is applied. 

In the development which follows, a perturbation analysis 
will be applied which when carried to the second-order will 
produce two linear boundary value problems from the one 
presently posed. Further, the boundary conditions of both 
linear problems will be referred to fixed boundaries. In 
addition, appropriate equations for determining force, and 
wave profiles will be obtained. 



i \ 
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B. PERTURBATION ANALYSIS 



First It Is assumed that $ and n each possess a perturba- 
tion series representation as follows 

( 10 ) 
( 11 ) 

where e is a small non-dimensional parameter which is a 
measure of the size of the ship's motion and is defined as 
e = a/b ("a" is the amplitude of the ship's motion and "b" 
is the ship's half-beam). Also in equation 10, and 

are defined to be the first and second-order velocity 
potentials, respectively, and similarly for and 

Note that as e-^0 there is no disturbance and therefore there 
are no zero-order functions. Substitution of equation 10 
into equation 2 yields 

+ O(e^) = 0 . (12) 

Equation 12 Implies that 



= 0 , and 


(13) 


= 0 , in R . 


(li<) 



A similar result holds for equation 9. 



$ = + O(e^) , 



n = en 



+ O(e^) , 
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1. The Linearized Free Surface Condition 



The following development of the linearized free 
surface conditions follows in general that of Stoker [353. 

The final linearized equations for the free surface boundary 
condition are obtained in the following steps. Equations 5, 
10, and 11 are substituted into equation 3, noting that p = 0 
on the free surface. Also, it is observed that 
and may be represented by the following form of the 

y 

Taylor series 



$^(x,y,t) 



= $^(x,0,t) 



2 vyy 



+ n(x,t)$^y(x,0,t) 
(x,0,t) + O(n^) 



( 15 )^ 



After collecting coefficients of like powers of e in the 
resulting equation, the following set of equations results 



gn 



( 1 ) 




0 



(16) 



gn 



( 2 ) 



+ $ 



( 2 ) 












on S^ . (17) 



Applying a very similar procedure to equation A yields the 
following result to the second order 



^The subscript v denotes partial differentiation with 



respect to x. 



or t . 



\ 
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1 



0 



( 18 ) 



n 



(1) - $(i) 

t y 



and , 



n 



( 2 ) 

t 






( 2 ) 



X X 



yy 



0 , on . (19) 



Differentiating equation l6 with respect to time, and 
eliminating from equations l6 and l8 yields the well 
known linearized free surface condition 

g$y^^ + = 0 , on y = 0 . (20) 



In equation 20, because of the use of the Taylor series 
expansions in n, the functions are evaluated for y = 0, 
rather than on the actual free surface. 

Following the same procedure, as before, for equations 
17 and 19 yields 




+ $ 



(2) ^ 1 

tt g 




($ 



( 1 ) 

tty 



+ 






( 1 ) 

yy 



) 



- 2(<J> 



X xt 



+ <I> 



y yt ^ 



( 21 ) 



Note the nonlinear function of on the right-hand side of 

equation 21. 

2 . The Ship Interface Condition 

Since the ship boundary is also in "motion'', a 
similar procedure must be applied there to "fix" the boundary 
S 2 . The development follows to some extent that of Potash 
[ 31 ]. Consider figure 3 in which the general ship hull form 
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FIGURE 3. Hull Coordinate System 

is described in two-dimensional cartesian coordinates in 
terms of the parameter s (arc length) in the follov;ing way 

X = x(s) , y = y(s) . (22) 

The inertial coordinate system is the same as shown in 
figures 1 and 2. The parametrically represented body 
coordinates (equation 22) correspond to the inertial coordi- 
nates of the body when the ship is at rest (neutrally buoyant 
position) . The motion of the ship in this reference frame 
can be defined in the following way 



x(s ,t ) 


= x(s) , 


(23) 


y (s ,t ) 


= y(s) + e y^^(t) , 


(24) 
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where the function is defined to be 

= Re{be^^^} . (25) 

In equation 25, a is the circular frequency of the harmonic 
motion and i is the imaginary unit. 

The velocity vector for the ship referred to the 
inertial reference frame is 

V = e yj^ J . (26) 

The dot denotes ordinary time differentiation. The unit 
outward normal vector (from the region) in terms of the 
parametric coordinates becomes 

n = -y’l + x'j , ( 27 ) 

where the prime indicates differentiation with respect to 
the variable s. Substituting equations 26 and 27 into 
equation 7 yields the ship-fluid Interface boundary condition 

^4>*n = eyj^x' . (28) 

Now expanding $ in a Taylor series in e (similar to free 
surface expansion) yields 

4>[x(s ,t) ,y (s ,t) , t ] = 4>(x,y,t) + (x ,y ,t ) + O(e^) . (29) 

Using equations 10 and 29, the normal derivative of $ on S 2 
may be expressed as 



' \ 
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y 



X y 

+ + x'<i>^^^ - y'yh<i>^y^ + x'y^$^y^]e^ + oCe^) 

on Sg . ( 30 ) 

Substituting equation 30 into equation 28 and equating like 
powers of e yields the following set of equations valid on 
the fixed boundary S 2 : 

*n = y^x* , 

^ y’5'h*xy\ 

Equation 32 may also be written 

= -y^(ii-V<J>^^^ )y . (33) 

It remains to consider the boundary (radiation 
boundary) and the function For consistency, must 

also be expanded in a perturbation series 

+ O(e^) . (34) 

Potash [31] and others have shown that and for 

a region of infinite depth, each represent a simple harmonic 
traveling wave of appropriate frequency and wave length. 
Zienkiewicz and Newton [43] have provided the appropriate 
homogeneous boundary condition for the surface Sj^ to success- 
fully pass a wave form without reflection. The equations are 



(31) 

(32) 
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.n = - -i- 

t 


y 


(35) 




C 2 t 




(36) 


In the 
defined 


above equations c^ and C 2 
by the following Implicit 


(the wave celerities) 
relationships 


are 




C-, = (g/o) tanh (^) 

Cl 




(37) 




Cp = (g/2a)tanh(^^) 

C 2 


• 


(38) 


3. 


Perturbation Formulas for 


Pressure , 






Force and Wave Amplitude 







The determination of dynamic loads on the ship 
requires the following equations which follow from the pertur- 
bation series assumed for $ and n 

p(x,y,t) = p^°^(x,y) + ep^^^(x,y,t) + e^p^^^ (x,y ,t) + O(e^) , 

(39) 

F(t) = eP^^^(t) + e^F^^^(t) + O(e^) . (40) 

In equation 40, F(t) Is defined as the total force per unit 
length required to sustain simple harmonic vertical trans- 
lation (see figure 3). 

Following a similar procedure as before, evaluating 
the Bernoulli equation (equation 3) on the ship's hull using 
equations 23, 24, and 29 yields 



' \ 
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- ^ p[x(s) ,y (s) ,t ] = gy(s) + e[gy^(t) + (x,y ,t ) ] 




By comparing equation 39 with equation we find 



p^°^ = - pgy(s) 



(^ 2 ) 



P^^^ = - PCgy^^ + y 



(^ 3 ) 







At this point observe that the second-order pressure is made 
up of contributions from both the first and second-order 
potential solutions. 



the force per unit length acting on the cylinder. However, 
in this treatment, it is desired to obtain relations for the 
exciting force as previously defined. The motive is that the 
theoretical result obtained should correspond to that 
measured by an appropriate physical experiment. The relation 
between the exciting force F, harmonic displacement yj^ and 
hydrodynamic pressure p, is 



Potash [ 31 ] and others have provided relations for 



F(t) = ep2Ay^ - / pn*j^ ds 



(^ 5 ) 
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where il Is one-half of the total symmetric arc-length of the 
ship’s hull (see figure 3). It should be noted here that 
Potash [ 31 ] has examined the question of the time-dependent 
total arc-length of the actual wetted hull and its effect on 
the limits of integration in equation ^15. However, it will 
be stated without proof that, provided the sides of the ship 
are vertical at y = 0, the second-order effect is zero. The 
only hull form examined to the second order satisfies this 
stated requirement. Therefore, equation 45 is exact for the 
case considered. 



considered since it is not in the definition of P. Accordingly, 
the following expression for P(t) results using equations 27, 

39, 43 , 44, and 45 



Comparing equation 46 with equation 40 yields the following 
set of equations 



In the calculation of the force F, the constant 
hydrostatic lifting force which comes from p^*^^ will not be 





(46) 



p(l) 




(47) 




(t) = p / + <I>^^^]x’ ds . 

(48) 
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V 



The wave amplitude may be written by inspection of 



equations l6 and 17 . 



n 



( 1 ) 



(x,t) 



1 

g t 



(^ 9 ) 



ri^^^(x,t) = 



- 

S t 



ty 









y 



(50) 



The boundary value problems for and may now 

be formally stated. However, at this point, a standard 
separation of variables technique is applied to eliminate 
time from the solution (since the solution is assumed 
periodic) . The appropriate definitions are 



(x,y,t) 


TD /a(I)/' \ icrti 

= Re{({)" (x,y) e } , 


(51) 


(x,y ,t) 


= Re{ (x,y ) . 


(52) 



Note the newly defined spatial potential functions and 

( 2 ) 

<t> are in general complex. Also because of the definitions 
in equations 51 and 52 it is necessary to replace y^ by the 
complex form 

yj = be^'^'*^ . (53) 



' \ 
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C. THE BOUNDARY VALUE PROBLEM FOR 

The boundary value for is formally stated using 

equations 9, 13, 20, 25, 31, and 35. 



= 0 , in R , 



- a2<j,(l) = 0 , on 



^(|) 



^ -n = iobx * , on 



^^(1) ->■ ia<j)^^^ 

-n = — ^ , on 



f 1') -»■ 

V4>' '.n = 0 , on and 



( 5 ^) 

(55) 

(56) 

(57) 

(58) 



D. THE BOUNDARY VALUE PROBLEM FOR 

( 2 ) 

Similarly, the boundary value problem for <t> may be 
formally stated using equations 9, 1^, 21, 25, 32, and 36. 



y2({,(2) = 0 , in R , 



(59) 



g4> 



( 2 ) 



( 1 ) 



. ,o2^(2) ^ (1) , (1)^ 



2g 



yy 



- , on S. 



^4) 



( 2 ) 

•n 



by'(j)^^^ bx'4)^^^ 

^xy ^yy 



, on S 2 



( 60 ) 



( 61 ) 



The factor h on the right-hand side of equations 60 and 
61 is due to complex algebra involved and is explained in 
Appendix A. 
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^<j) 



( 2 ) 



n = - 



21a({) 



( 2 ) 



on Si 



( 62 ) 



= 0 , on and Sj 



( 63 ) 



Note the striking similarities between the first and second- 
order problems. Note also, that the nonlinearities in the 
original problem have been approximated by the nonlinear 
terms in in equation 60, and the nonlinear terms in 

and 4)^^^ in equation 6l (y^^ does not explicitly appear). 

2 

The additional nonlinear amplitude dependence comes from e 
in the original perturbation expansions (equations 10 and 
11). The structure of the problem requires the solution of 
first in order to define the solution for 
The conversion • to complex algebra allov;s another observa- 
tion to be made. Since only the real part of is 

meaningful. Appendix A demonstrates the required complex 

manipulation. Close examination of that result indicates 

( 2 ) 

that the nonhomogeneous -boundary conditions on <t> contain 

( 2 ) 

time-independent terms. This implies 4> has a time- 
dependent and a time-independent part. However, since the 

physical quantities of interest in the solution all involve 

( 2 ) ( 2 ) 
derivatives of <{) , the time-independent solution of 4 

will not be pursued further. The above comments are also 

true for the second-order functions of pressure, force, and 



5 

terms 



Lee [l8] comments on the time-independent solution in 
of mass transport. 
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wave amplitude. However, the time-independent quantities 
are of interest and have physical meaning as will be seen 
later. 

( 2 ) 

The boundary value problem for (j) is valid only in 
\ 

water of infinite depth. The following discussion will 

explain the problem and develop a boundary value problem 

( 2 ) 

which will allow for the solution of (j> in finite depth. 

E. THE SECOND-ORDER BOUNDARY VALUE PROBLEM 
WITH FINITE DEPTH 

Potash [31] has shown that the asymptotic form of the 
velocity potential $ as x->o° becomes (in water of infinite 
depth) 

(610 

Where denotes the wave number of the first order wave 

A 

(k^ = a/c^) and and G2 are coefficients which depend on • 
frequency a. However in the more general case of finite 
depth, an additional term due to the Stokes'^ second-order 
potential appears as follows 

./ n o / u i[at-(knx+Y^^^)] . 2 r„ -l(koX+Y^^^) 

= ReleH^e +e [H2e ^ d ' 

- i H^^cosh(2kih) -2i (kix+Y ^ ^^)-i 2iat, . (65) 

8 c^^ slnh^(k^h) ® J e J 



'^The Stokes potential to the second order is v;ell known 
to be zero in infinite depth. 
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In equation 65, ^2 is defined as the second-order wave number 
(k 2 = 20 / 02 ). ^1 ^2 constants similar to and G 2 . 

Close examination of equation 65 reveals that the radiation 
boundary condition (equation 62) would fall since the 
"Stokes wave" which Is present Is traveling at the same 
celerity as the first-order solution. This means that the 
radiation boundary would "reflect" the Stokes’ wave since 
the boundary condition defined by equation 62 will only pass 
a wave traveling at a celerity of C 2 , defined by equation 38 . 
A successful resolution of the above described difficulty 
requires the formulation of a different boundary value 

(2) 

problem. We define the complex velocity potential ip as 

^.( 2 ) - 13 a _ (65^) 

8c^^ slnh^(k^h) 



which corresponds to the third term In equation 65 . Then, 

following previous conventions, the complex velocity potential 
—( 2 ) 

(|)^ Is defined as 






( 66 ) 



1. The Second-Order Boundary Value 

ZTpl 

Problem for ' 

(o\ 

The definition of ' In equation 66 allows the 

formal statement of the following boundary value problem, 

( 2 ) 

which follows from equation 66, the definition of \p and 

( 2 ) 

the boundary value problem for 4> • 



il2 



I 






In R 



(67) 



V2 ^(2) = 0 



,-(2) „„2-(2) _ , r ^2,(1) ^ ^,(1) 



^ - 4a ^ = i 






X y 



gt|i^2) ^ ^ on S^, 



( 68 ) 






b<}) 



( 1 ) 



2) = yt ( ^ L- + ,|j(2)) 



b<}) 



( 1 ) 



. P(:zz^. ,(2)) , 



2 ■ ■*'y ^ > on S^, 

(69) 






21 a(j) 



( 2 ) 



, on S, 



(70) 



;^-(2) -> 
V(})'' *n = 



-*■ (2) ->■ 

- Vi|j^ ' *n , on S. 



(71) 



‘11 = 0 , on S 



5 • 



(72) 



Equation 70 will now properly behave because the asymptotic 
— ( 2 ) 

form of <}) using equations 65 and 66 Is 

*(2)(«,y) e21°'= = Hj ^ 2 ±ot (^3, 

/ 2 N 2 i 0 ^ t 

Note the celerity of the complex wave form ' e‘^ Is C 2 

and the radiation boundary condition will now work. 

F. FORMULAS FOR PRESSURE, FORCE, AND WAYE AMPLITUDE 

The development which follows requires additional defini- 
tions, necessary because of the separation of variables scheme 
previously Introduced. They are: 



^3 



i 



F^^^(t) = Re{f^^^e^^^} , F^^^(t) = Re {f } , 

n^^^(x,t) = Re{n^^\x)e^°^} , = Re{n^^^x)e^^^^}. 



The functions p^^\ P^^\ are complex. This 

also true for the constants f^^^ and f^^\ 

Using the above equations and equations 43, 44, 47, 
49, and 50 the equations for pressure, force, and wave 
height for both first and second-order solutions follow 



= -p[bg + io4)^^^ ] , 



p(2) = 






.( 1 ) _ 



p ^ ( 1 ) 

■2pbAa + p / (gb + ic(j)^ ')x' ds , 



~Z 



( 1 ) 



(2) = .p . 21o^<2>) 



iT'^> = 



^(2) = 



la4) 



( 1 ) 



g 



1 ( 2 ) 

- - f2ia(f)^"^'' + ^ ^ — 

g ^ 5 2 2 2 



]) 



(74) 

(75) 

(76) 
is 
48, 

(77) 

(78) 

(79) 

' ds , 

( 80 ) 

( 81 ) 

( 82 ) 



Again, equations 78, 80, and 82 include the special considera- 
tion of the complex algebra involved which is explained in 
Appendix A. 



III. FINITE ELEMENT DISCRETIZATION 



The Finite Element Method and its application to scalar 
field problems is described in various texts, including one 
by Zienkiewicz [^5]. The essential ideas involved follow. 
The region of figure 2 is replaced by a grid or mesh of 
finite elements as shown in figure ^ . 



FIGURE 4. Finite Element Mesh 




Then it is assumed that the field function sought can be 
represented by the following equation 

<l)(x,y) = N ^ , in R . (83) 

The total number of nodes in the region will be defined to 
be m. Then N is a 1 x m row vector of shape (interpolating) 
functions and <{> is an m x 1 column vector of nodal values of 
the complex potential <}) (in this problem) . Equation 83 
implies that knowledge of the vector <p defines 4 >(x,y) 
everywhere in R. 

I V 
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There are several choices of shape functions which may 
be used In equation 83 . This work selected the cubic shape 
functions described by Zlenklewlcz [^5] as the "serendipity 
family". The elements used are Isoparametric. A brief 
description of the 12-noded Isoparametric elements Is given 
In Appendix B for the unfamiliar reader. The key advantage 
of the Isoparametric element lies In the fact that regions 
with curved boundaries may be readily represented. 

The development that follows defines the set of linear 
equations which determine the vector (j) of equation 83 . 

One general approach to discretization of several types 

of field problems In finite element analysis Is to define a 

functional whose Integral over the domain of the problem Is 

to be minimized. The Calculus of Variations shows that If the 

Euler equation of the Integral to be minimized Is the governing 

equation (l.e., Laplace's Equation) then the two problems-are 
7 

equivalent . 

However Zlenklewlcz and Newton [43] have shown the 
discretization may be readily accomplished by applying the 
Galerkln process. The weight functions are chosen to be the 
nodal shape functions. Mathematically stated 

/ + (j) ) dR = 0 . (84) 

XX J J ^ 



7 . . 

This statement assumes the same boundary conditions are 

used. 



The superscript T denotes transposition and the function (j) 
is considered to be any of the complex potential functions 
whose solution is sought. Application of the divergence 
theorem transforms equation 8^ into 



! [N^ 

R 





dR 



/ N V(j)*n dS , 
S " 



where S denotes the boundary of the region R. 



( 85 ) 



A. THE DISCRETIZED FIRST-ORDER PROBLEM 

Application of equations 8^1 and 85 to the first-order 
boundary value problem stated in Chapter II yields 
the following matrix equation 



(-a^Q + — D + H) = labr^^^ . (86) 

. 1 ^ 

The definitions of the matrices in equation 86 follow. The 
matrix H is m x m, real, symmetric, and defined as 



H = 



/ 

R 






dR . 



( 87 ) 



The matrix H comes from the Integral on the left-hand side 
of equation 85 . The matrix D is mxm, real, symmetric, and 
defined by 



D 



/ 

S 



n'^N dS 



( 88 ) 



^7 
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D is contributed by the boundary integral on the right-hand 
side of equation 85 and represents the discretization of the 
radiation boundary condition (equation 57). 

The matrix is real, -m x m, symmetric and given by 

= ^ / n'^N dS . (89) 

S3 

parallels D in its origin and represents the contribution 
from the homogeneous free surface boundary condition (equation 
55) . 

The vector r^^^ is real and defined by 

= / n’^x' dS . (90) 

The above vector is the discretized ship-fluid interface 
kinematic boundary condition (equation 56 ). The surfaces 
S 3 , S^j, and have homogeneous boundary conditions and 
therefore make no contribution to r^^^. It follows that r^^^ 
is non-zero only at nodes along the ship's hull (surface S 2 ) . 

B. THE DISCRETIZED SECOND-ORDER PROBLEM 

—( 2 ) 

The problem for only will be presented since it 

( 2 ) 

contains the solution for (j) in the limit. The boundary 

— ( 2 ) 

value problem for <}) yields a linear system of complex 
equations similar to equation 85 . 






2io 



D + H ) (p 



( 2 ) 



( 2 ) 

_ 






i 



(91) 



In equation 80 receives 



( 2 ) 

The right-hand side vector r 
contributions from the surfaces S^, S2 and S^. Accordingly 
a further decomposition Is necessary; 





+ 




+ 




( 92 ) 



( 2 ) 

The vectors r^ are defined below. First to be consistent, 

( 2 ) —( 2 ) 
the function \p must be represented In the same way as c]) 



^ ^ ^( 2 ) 



( 93 ) 



(2) 

Then r^ becomes, from equations 71 and 85, 



( 9 ^) 



r^^^^ = / n'^N dS . 

^1 

(2) 

The nodal values of if; are determined by partial dlfferen- 

'S' X 

( 2 ) 

tlatlon of the function if; and then representing the 
( 2 ) 

function if; In similar manner as equation 93 . 

X 

( 2 ) 

The vector r^ ' Is defined by equations 69 and 85 



2 ix 



<f>(l)b 



— + if;^^^]dS - f — + if;^^^]dS 



( 95 ) 



( 2 ) 

similarly r^ ' follows from equations 68 and 85 



r(2) = 1 , N^H[- I” 

-.3 Ssj — 



X Q 

2 g ~ 2 ^ 



lOY - + ^a^if;^^^] dS . 

^ ^ u ^ 

( 96 ) 
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The 1^ component of a Is 






( 97 ) 



where and are the 1^^ components of and 

1 ijy ~ ~y 

respectively. Similarly 






1 i,yy 



and 



( 98 ) 



Y, = 









( 99 ) 



The nodal values of the various partial derivatives of 
are determined by differentiating equation 83 . The method 
of calculation is discussed in the next section. The reader, 
unfamiliar with isoparametric elements, should refer to 
Appendix B prior to reading the next section. - - 

1 . Calculation of Nodal Values of the 
Field Derivatives of 4> 

The vectors ^vv^ required in 

equations 95 and 96 are determined by applying standard 
concepts of finite element analysis. The appropriate relation 
for the l'^ node of a given element is 



40 e 


= j ^ 
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1 

0) 

1 



(l)e 



( 100 ) 



50 



l \ 



J Is the 2x2, real Jacobian matrix of the transformation. 

The vector N Is a 1 x 12 row vector of element level shape 

functions. The superscript "e" denotes all quantities In 

brackets are at the element level. In equation 100, the 
-1 e 

elements of J and N are evaluated at node 1. 

Application of equation 100 to all twelve nodes of a 

(l)e (l)e 

given element determines the vectors ((> and (J)„ 

(l)e (l)e 

Then, replacing (|) by (j)y In equation 100 yields the 

CDe d')e 

components of the vectors <|) and . Mathematically 

'^'Xy *^yy 

stated 







nJ 


^i,xy 


= r" 








N® 






~n 



. ( 101 ) 

^ *y 
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IV. THE COMPUTER PROGRAM AND NUMERICAL SOLUTION 



A. MAIN FEATURES OF THE COMPUTER PROGRAM 

The formation and solution of the linear system of equations 
defined by equation 86 ((|>^^^ system) or equation 91 
system) was accomplished on an IBM 36 O /67 computer at the Naval 
Postgraduate School. The computer has a maximum word length 
of sixteen bytes. 

A comparison of equations 86 and 91 demonstrates the close 
similarity of the coefficient matrices for the first and 
second-order solutions. Note that only the coefficients of 
Qq and D change. An additional observation follows from 
examination of equations 87 , 88, and 89 defining H, D, and 
respectively. Note that all three matrices are purely a 
function of mesh geometry. Therefore, once a mesh is chosen 
to define the region of interest, the three matrices previously 
mentioned need only be calculated once regardless of the forced 
motion . 

A well known property of the finite element method is that 
it usually produces a banded, symmetric, coefficient matrix. 

That is the case here and this has obvious computational 
advantages. However, in this problem, the final system of 
equations is a function of the frequency of motion as for 
example in the first-order problem 






iabr 



( 1 ) 



( 102 ) 
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The matrix Is an mxm complex matrix defined by equation 

86 to be 

= -a^Q^ + — D + H . (103) 

- ^ ^ ^ 

Although is complex and frequency dependent, it is 

possible to perform most of the triangular decomposition of 

(by Gauss elimination) only once, using real arithmetic, 
for a given problem geometry. This economization was first 
employed by Chenault [4] and is described in more detail in 
Appendix C. In addition, this author employed a vector 
storage of the upper triangular portion of which greatly 

reduced the storage requirements of the program. The vector 
storage scheme made it possible to store only the portions of 
which are non-zero. 

g 

The net advantage of employing Chenault *s scheme and the - • 
compact ' vector storage was that the system of 733 complex-, 
equations which is generated by the mesh shown in figure 4 
could be solved for both the first and second-order problem 
twenty times in twenty minutes. Further, only ten percent of 
the total matrix was actually stored and the computer solution 
was accomplished in core. 



O 

°The Chenault scheme described in Appendix C was developed 
by Professor R. E. Newton at the Naval Postgraduate School. 
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B. NUMERICAL SOLUTIONS 



The solution of the first-order problem for a given mesh 
geometry is required before the second-order solution can be 
obtained. Further, the first-order solution must be highly 
accurate if a good second-order solution is to be expected. 

This fact is made more obvious by re-examining equations 68 
and 69 in part II. 

1 . First Order Numerical Solutions 

Several first-order problems were solved and closely 
examined for accuracy. The numerical solutions are presented 
in part V. 

a. Mesh Requirements 

One of the primary considerations in choosing a 
proper mesh concerns the positioning of the radiation 
boundary (S^^). It may be recalled that the surface S^ must 
be far enough away from the ship hull (location of disturbance) 
so that the local effect of the disturbance has decayed and 
only an outgoing wave is present. Bal [3] has addressed this 
question analytically. He has shown that the eigenfunctions 
associated with the non-propagating velocity field decay 
exponentially with x. Bai recommends the criterion 

w = ah + b , (10^) 

min * 

where a is a monotonically increasing function of frequency. 

The range of a is from 1.5 (low frequency) to 3- The values 
of a were determined based on the decay rate of the most 
persistent eigenfunction and an attenuation factor of 0.01. 
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This investigation has found the above criterion somewhat 
conservative. A value of a = 2/3 was found to be adequate 
for the velocity potential to numerically approach its 
asymptotic limit for all frequencies studied for the semi- 
circular hull (h/b = 6). Similarly, a value of a = 9/10 was 
found to be adequate for the bulb hull form described in 
part V (h/b = 10) for all frequencies studied. The values 
of a specified above were determined for the deepest regions 
studied for each hull form. The region widths used in this 
case (w/b = 6, w/b = 10) were found to be more than adequate 
for shallower depths. 

Another requirement on w which tends to restrict 
its maximum possible value comes from the practical limits on 
mesh fineness and the resulting computer storage requirements. 
The mesh must be sufficiently fine to represent properly a 
traveling surface wave. The finer the mesh, the larger the 
computer storage requirements and the longer the computer 
solution time. Therefore, this latter restriction tends to 
oppose the former one. 

Bal [ 3 ] used quadratic elements and recommended 
at least 10 surface nodes per wave length. Visser and Van 
der Wilt [39] also used quadratic elements and recommended 
that an element could not span more than one quarter of a 
wave length. This investigation examined this question in 
detail. As previously mentioned, cubic isoparametric elements 
were used. It was found that if the gravity wave did not have 
to be propagated more than two wave lengths, the elements 
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could be allowed to span six-tenths of a wave length. 

However, as in the case Visser and Van der Wilt [39] studied, 
the successful propagation of a gravity wave over many wave 
lengths requires a finer mesh spacing. The one quarter wave 
length Visser and Van der Wilt recommend was found to be a 
practical upper limit. 

A parallel problem to positioning the radiation 
boundary concerns the location of the bottom (surface S^) 
to simulate infinite depth. Comparison of the finite element 
solutions obtained in this study with analytic solutions 
based on infinite depth led to the conclusion that a keel to 
bottom clearance (h - d) of 5b will adequately simulate 
infinite depth. However, satisfying the classical hydro- 
dynamic requirement that the depth must be at least one-half 
of a wave length becomes impractical at low frequencies. 
Ignoring this requirement did not seem to affect solution 
accuracy . 

b. Mesh Data Preparation 

The meshes for a given boundary geometry were 
prepared using a very versatile mesh generator computer 
program developed by Adamek [1] at the Naval Postgraduate 
School. This tool greatly expedited mesh preparation. For 
example , the entire data for the mesh shown in figure ^ 

(733 nodes, 132 elements) could be generated in twenty 
seconds. A modification was used to supply accurate values 
of the nodal coordinates along the ship-fluid interface . 
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c. Convergence of First-Order Solutions 



Convergence of the first-order solutions could 
not be determined by standard numerical techniques (i.e., 
mesh refinement) due to computer size limitations. However, 
the excellent agreement between the first-order solutions 
obtained in this study and both theory and experiment left no 
doubt that the solutions were quite good. 

2 . Second-Order Numerical Solutions 

The successful solution of the first-order problem 
made study of the second-order problem possible, 
a. Mesh Requirements 

All of the requirements on mesh geometry previously 
mentioned for the first-order solutions of course apply to the 
second-order solution. However, the second-order problem 
demands even more stringent requirements. One. reason for this 
follows from the fact that the wave lengths of the propagating 
portion of the second-order velocity potential are .much 
shorter (one-fourth of the first-order wave lengths in Infinite 
depth solutions). Therefore, since it is impractical to use 
a different mesh for each solution, the mesh must have a free 
surface element spacing approximately one-fourth of that 
required for the solution of for a given frequency range 

to be studied. The mesh shown in figure ^ was used to solve 
the Infinite depth second-order problem. It contains 132 
elements and 733 nodes. 
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b. The Second-Order Boundary Condition on 

Additional constraints on mesh geometry were 
found to be required. These constraints arise from the 
requirements to calculate second-order partial derivatives of 
which appear in equations 95 and 98. It was essential 
to have rectangular elements with equal node spacing along 
the free surface for success in representing the non- 
homogeneous boundary condition given by equation 96. 

The further discussion of this point will be 
aided by the definitions of some complex functions. Define 



Q(x,0) 

R(x,0) 



2g ^y 






(105) 

(106) 



In terms of these definitions the right-hand side of equation 
68 is 

P(x,0) = Q(x,0) - R(x,0) . (107) 

The reader should recall that the successful solution of the 

( 2 ) 

boundary value problem for tj) using the method of finite 

elements required the definition of a new boundary value 

—( 2 ) 

problem for the function (p • The reason was that in the 
finite depth problem a second-order Stokes’ wave is produced 
with celerity c^ as well as a second-order wave with celerity 
C .2 (contained in '). The presence of the two waves with 
different celerities meant the radiation boundary condition 



would fail. The problem was mathematically eliminated by 

subtracting the well-known analytic solution for the Stokes 

( 2 ) 

wave portion of the second-order solution of (|) on every 
boundary of the region R. The mathematical steps are straight- 
forward as shown In Chapter II. The asymptotic value of the 
function P(x,0) Is zero as x gets large. But, the function 
R(x,0) which represents the classic Stokes second-order wave 
free surface boundary condition has the form 

R(x,0) = — ^ U^^^(w,0)|^ g-2i(kix+Y^^^) ^ 

g'^sinh'^k^d 

everywhere on S^. Therefore, if the numerical representation 
of the function P(x,0) is to approach zero near the radiation 
boundary; the numerical representation of the function Q(x,0) 
must approach R(x,0) or 

Q(x,0) -»■ R(x,0) as x “ (109) 

Note that the amplitude and phase of must be known to 

( 2 ) 

completely define R(x,0) and therefore \p . 

The above discussion points out the necessity of 

a very accurate numerical determination of the function Q(x,0). 

— (2 ) 

The second-order solutions for (j) were closely examined near 
the radiation boundary and it was found that on the average 
the function Q(x,0) was within one-percent of the analytic 
function R(x,0). This accuracy was achieved even though 
formation of Q(x,0) requires calculation of second derivatives 

' \ 
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of using element shape functions which guarantee only 

C° continuity across Inter-element boundaries. 

c. The Ship-Fluid Interface Second-Order 
Boundary Condition 

The accurate representation of the ship-fluid 
Interface (S2) boundary condition for the second-order 
boundary value problem (equation 69) Is of course essential. 
Potash [31] encountered difficulties In making a numerical 
calculation on this boundary because the potential singulari- 
ties which define his solution lie along S2 . This results 
In the required functions and being only piecewise 

continuous. This problem required an additional approximation 
to be made. 



The method used herein also has the problem that 

x(l) and are not continuous across element 

» ^y ’ ^xy ’ ^yy 

boundaries. However, no singularities are present. Extensive 

Investigations In this work have shown that It Is essential 

that the elements along the hull be as square as possible, 

and further that the node spacing on the elements along the 

hull should be uniform. It was found that If the hull 

elements were distorted (not square or rectangular) the 

numerical representation of the functions and was 

highly unstable. However, the mesh requirements for 

numerically calculating and along the hull were 

X y 

much less stringent. 

The adequacy of the hull elements used In this 
Investigation (see figure was tested by defining a known 
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potential function on the elements along S2 



+ i xy 



(110) 



and then applying the computer program's numerical scheme 



It was found that the cubic, isoparametric elements along S2 
produced an average error of about one percent in calculating 



five percent at isolated nodes. 

d. Convergence of Second-Order Solutions 

Again, computer capacity limitations prevented 
the application of systematic mesh refinement to determine 
convergence of the second-order solutions. However, it was 
possible to refine the elements next to the ship's hull and 
no significant change in the solutions was observed (the 
elements along the free surface were performing properly) . 
Further, it is shown in Chapter V that second-order results 
for infinite depth are in satisfactory agreement with those 
of other researchers. No comparison is available for finite 
depth solutions, but a smooth transition is observed as the 
depth becomes effectively infinite (See Chapter V) . 





(Ill) 



and 




( 112 ) 
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V. NUMERICAL RESULTS 



A. DEFINITIONS 

The definitions which follow are useful for presentation 
of results. 

1 . Added Mass and Damping Coefficients 

The complex amplitude F^^^ of the first-order 
resultant hydrodynamic force acting vertically on the ship's 
hull is 



F^^^ = / p^^^x’ ds . 

^ -A 



( 113 ) 



The real part of F^^^ is in phase with ship’s acceleration 

O’ 

and the ship’s added mass is defined by 



^a = 



Re{F^^^} 



( 114 ) 



2 

where the ship’s acceleration amplitude is given by a a. 
The added mass coefficient is defined to be the ratio of 
added mass to ship’s displaced mass 



m Re{F^^h 

2pAaa^'“ 



( 115 ) 



The imaginary part of F 
velocity of the ship and is the 



is in phase with the 

y 

n 

”non-conservatlve"^ 



^The potential theory is, of course, conservative. However, 
the waves which propagate to "infinity" represent work done by 
the ship which is effectively lost. 
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component. Accordingly, the damping coefficient is defined 
to be 



C , = ^ 

2pAaa"^ 



(116) 



and are frequency dependent for a given hull form and 
region depth. 

A dimensionless frequency parameter is defined as 

follows 




The parameter 6 is commonly used by other researchers and is 
used herein for ease of comparison of solutions. The 
parameter 6 is equivalent to in infinite depth (X is the 

wave length) . 

The added mass and damping coefficients previously 
defined are for heave. The same coefficients for sway motion 
are similarly defined. The only difference in the definitions 
is is replaced by the horizontal force amplitude F^^^ 

y ^ 

defined by 



F 



( 1 ) = 
X 



f p^^^y' ds , 



(118) 



where p^^^ is the hydrodynamic pressure due to the swaying 
motion. There is a coupling between sway and rolling motion 
for a general ship cross-section. However, the only hull 
form studied in the sway mode was semi-circular in 
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cross-section and the coupling force for this case is zero. 
Therefore, the appropriate coefficients typically used are 
not defined herein. 

One further comment is warranted here. The second- 

order time-dependent hydrodynamic forces are harmonic with 

frequency 2a. Therefore the definition of second-order added 

mass and damping coefficients is meaningless . 

2 . Dimensionless Force Coefficients 

The force coefficients f^^^ and f^^^ are also 

presented in dimensionless form. The non-dlmenslonallzlng 

factor for the various force amplitudes is 1/pgb . Table I 

presents the correspondence between the dimensional and non- 

dimensional force quantities. The second-order quantities 
( 2 ) ( 2 ) 

fg and f^ which appear in table I are amplitudes of the 

time-independent (static) and time-dependent (dynamic) 

( 2 ) ( 2 ) 
portions of f respectively. The existence of f is .. 

s 

readily verified by examining equations 75 and 80 in light of 
the discussion in Appendix A. 



TABLE I 

DIMENSIONLESS FORCE COEFFICIENTS 

Force Amplitude f^^^ f^^^ f^^^ 

Dimensionless ^d^^ 

Force Amplitude 



3 . Wave Amplitude at Infinity 



Another suitable measure of the damping of ship's 
motion (to the first-order) is the wave amplitude at 
"infinity" (l.e., near the radiation boundary). Further, 
experimentally it has been found to be an easier quantity 
to measure (particularly for high frequency motion since the 
damping force is very small then). Accordingly, the dimen- 
sionless first-order wave amplitude at infinity is defined 
by 





x=w 



a 



(119) 



B . HULL FORMS 

Other investigators have studied a variety of hull forms 
both by experiment and analytically. The semi-circular hull 
wa's'“chosen because of the wealth of data available for 
comparison of both first and second-order solutions. The 
bulb-shaped hull studied by Paulling and Porter [29] was 
also examined for first-order solutions. The shape of the 
bulb hull in comparison to the semi-circular hull is shown 
in figure 5 • 

The coordinates of the bulb hull studied are defined by 
the following mapping 



X + ly 




( 120 ) 



\ 
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FIGURE 5 . Hull Forms 

The complex number z, represents the points on a circle In 
the ^-plane . The choice of radius In the ^-plane defines 
the scale of the hull cross-section. The values of the e^ 
are given In table II. 



TABLE II 

CONSTANTS FOR HULL FORMS 



Hull 


®1 




"5 


A/bd 


Seml-Clrcular 


0.0 


0.0 


0.0 


0.785^ 


Bulb 


-0.7132 


-0.02096 


0.0605 


0.6950 



C. FIRST-ORDER NUMERICAL SOLUTIONS 
I . Seml-Clrcular Hull 

The semi-circular hull form has been studied by many 
researchers mentioned previously In part I. Specifically, 
the original solution by Ursell [ 38 ] for heaving motion has 
been confirmed analytically by Porter [30], Vugts [^0], and 
others. Ursell ’s solution has also been experimentally 
confirmed by Paulllng and Porter [29] and Vugts. The 
analytic results of Ursell for and are shown by the 
continuous curve In figure 6. The numerical results obtained 
by the finite element method (PEM) In this work are repre- 
sented by the circled points In figure 6. The FEM points 
shown were obtained using a mesh similar to that shown in 
figure , but using only 62 cubic elements and 359 nodes. 

The extension of the solutions for C_ and C, to 

m d 

finite depth Is easily accomplished by simply redefining the 
mesh geometry (decreasing the vertical dimension) . Kim [lA] „ 
and Yu and Ursell [42] reported analytic results. The .. c_.,, 
results obtained by the FEM in this work for h/d = 4 agree 
with those reported by Kim for and within three 

percent over the range 0.25 < 5 < 5- The results of Kim are 
in disagreement with those obtained by Yu and Ursell for 
h/d = 2. The results obtained herein fall between the two. 

Bai [ 3 ] reports inability to verify Kim's result, but gives 
no details. Bai used the finite element method to obtain 
his solutions as mentioned in part I. 
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FIGURE 6. Added Mass and Damping Coefficients in 

Heave, Semi-Circular Hull, Infinite Depth 

I 

Figure 7 shows the effect of finite depth on C^^ 
obtained in this work. The figure shows that significant 
changes do not occur until very shallow water is encountered 
(h/d _< 2) . Further, the effect of shallow water on C^^ is 
more pronounced at very low and very high values of 6 . Note 
the unusual variation with depth for 6 = 0.5. decreases 

with decreasing depth (opposite the trend at higher values 
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1.25 




FIGURE 7. Variation of Added Mass Coefficient with 
Depth, Semi-Circular Hull in Heave 



of 6). Although it does not appear true in figure 7, 

(for 6 = 0.5) has reached its asymptotic value at h/d = 6. 
Figure 8 shows the corresponding effect of finite depth on 
C^. The figure indicates that is only sensitive to 
shallow water effects at very low values of 6. 




h/d 



FIGURE 8. Variation of Damping Coefficient with 
Depth, Semi-Circular Hull in Heave 
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The semi-circular hull was also studied in sway 
motion to the first-order. The only changes in the boundary 
value problem for defined in part II are in equation 56 

which defines the surface . The changes are 

V({)^^^-n = labx' , on , (121) 

for equation 56, and 

=0 , on , (122) 

for equation 58 . The boundary condition on defined by 
equation 58 remains unchanged. The ship’s motion Is defined 
in a manner analogous to equations 23 and 24 

x(s,t) = x(s) + e > ( 123 ) 

y(s,t) = y(s) . (124) 

The function x^(t) Is defined the same as yj^(t) 

i^t 

= Re{be , (125) 

and e is defined in the case of sway motion to be 

e = I ' (126) 

The symbol "a" is the amplitude of the sway motion in this 
case . 

Vugts [40] presents both experlm.ental and analytic 

solutions for C and C, for the semi-circular hull in sway. 

m d 

Figure 9 shows a comparison of Vugts’ theoretical solution 
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FIGURE 9. Added Mass and Damping Coefficients in Sway, 
Semi-Circular Hull, Infinite Depth 
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with the FEM results obtained herein. The FEM calculated 
points are shown by the solid circles and squares In the 
figure. The solid curve (from Vugts’ theory) was scaled 
from a very small graph and may be somewhat In error. 

2 . Bulb Hull 

Paulling and Porter [29] reported both analytic and 
experimental results for the bulb hull form. They show 
excellent agreement by comparing first-order exciting force 
amplitude coefficients F^^^ and wave amplitude coefficients 
respectively. Figures 10 and 11 compare the values of 
and obtained by the FEM herein (circled points) 

and Paulling and Porter’s theoretical calculations (solid 
line). The FEM mesh parameters were h/d = 2 and w/b = 25. 

This value of w/b was found later to be quite conservative 
(w/b = 10 would have been adequate) . 

Figure 12 shows added mass and damping coefficients 
for the bulb hull vs h/d for various values of 6. The curves 
have a similar trend to those of the semi-circular hull form. 
The corresponding curves of force F^^^ and wave amplitude 
are not presented. However, the maximum change in 
exciting force amplitude does not exceed five percent for 
h/d = 1.2 as compared to infinite depth (h/d = 2). Similarly, 
the maximum change in was only ten percent over the same 

range of h/d. 

3 . Wave Generator Studies 

MacCamy [20], in 1957, reported an analytic theory 
using a Green's function approach for solving the problems 
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FIGURE 10. Exciting Force Amplitude in Heave, 
Bulb Hull, Infinite Depth 




FIGURE 11. Wave Amplitude for Heave, 
Bulb Hull, Infinite Depth 



0.4 




FIGURE 12. Variation of Added Mass and Damping 

’ Coefficients with Depth, Bulb Hull in Heave r 



of gravity wave production by a moving partition and gravity 
wave reflection from a horizontal strip. The problem of 
gravity wave production by a moving partition was studied in 
this work to further help determine proper mesh characteris- 
tics for gravity wave propagation. 

The wave generation problem is readily placed in the 
framework of the class of problems being studied herein with 





I 






m 




only minor modifications.^*^ All of the boundary conditions 
posed for in part II remain unchanged except those on 

and S2 . The surfaces and S2 become a vertical line. 
The region therefore is rectangular. The left-hand boundary 
will be designated as as shown in figure I3. 




The boundary condition on becomes 

= q(y) , on SJ^ , (127) 

The function q(y) describes the mode of displacement of the 
wave generator. The plunger type and flapper type (shown in 
figure 13) wave generators were studied and compared with 



This observation is not true in the potential theory 
singularity distribution approach. The source singularity 
must be replaced by doublets and the problem completely 
reformulated . 
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MacCamy's results. The function q„(y) for the plunger type 
generator Is simply 



where is the maximum amplitude of the wave generator 
motion. Similarly, the function q^^Cy) for the flapper type 
generator is 



The comparison of MacCamy's results with those obtained 
herein is shown in figure 1^ for the plunger. The circled 
points were obtained by the FEM and the solid line is MacCamy's 
theoretical curve. Similarly, the comparison of MacCamy's 
result for the flapper type wave generator vs the FEM results 
' obtained in this work is shown in figure 15. The excellent 
^agr-eement allows two observations to be made. First, the FEM ^ 
approach to solution of problems of the type being presented 
is highly versatile in its application. Second, the radiation 
boundary placement was further verified. 

D. SECOND-ORDER NUMERICAL RESULTS 

1 . Semi-Circular Hull — Infinite Depth 



confidence that second-order solutions could be attempted. 
The only hull form studied, to the second-order, was semi- 
circular. The accuracy of the second-order solutions was 



qpO) = lip 



( 128 ) 




(129) 



The first-order solutions just presented established 
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FIGURE 1^. Wave Amplitude Generated by a 
Plunger Type Wave Generator 
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FIGURE 15. Wave Amplitude Generated by a 
Flapper Type Wave Generator 
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established by comparing with the work of other researchers, 

Lee [ 17 , 18 ], Parlssls [27], and Potash [ 31 ] have reported 

numerical solutions for second-order exciting force coeffl- 
—( 2 ) —( 2 ) ( 2 ) 

dents , F^ and phase angle . Figure I 6 compares the 

—( 2 ) —( 2 ) 

values of F^ and F^ reported by Lee, Parlssls, and 
Potash with the numerical results obtained by the FEM In this 
work. The points presented from the other researchers' works 
were all obtained from Potash [31]. Potash Indicates Lee's 
points are corrected In some manner from his original work 
[ 17 ] and otherwise unpublished. 

f p ) 

Figure 16 shows excellent agreement for F^ ’ and very 

s 

2 ) 

good agreement for F^ . Only Potash's points show signi- 
ficant deviation and then only for 6 near 1.5. Potash 
reported some numerical difficulties for a range of 6 near 1 
and therefore Potash's results are not presented for those 
values of ' 6 which are suspected to be In error. Similarly, 
figure 17 shows a comparison of ' (phase angle) results. 
The figure shows the agreement Is not as good. However, none 
of the other works shown (Lee, Parlssls, Potash) agree In any 
consistent fashion over a large range of 6. The reason for 
the discrepancy has not been resolved. 

2 . Seml-Clrcular Hull — Finite Depth 

— f 2 ) 

Figure 18 presents the effect of depth on F^ and 

as obtained by the FEM In' this work. The author knows 
s 

of no other solutions. Solutions were obtained for values 
of h/d from 6.0 (Infinite depth) to 1.5. The figure shows 
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FIGURE l6 . Second-Order Exciting Force Amplitude in 
Heave, Semi-Circular Hull, Infinite Depth 
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FIGURE 17. Second-Order Exciting Force Phase Angle 

in Heave, Semi-Circular Hull, Infinite Depth 
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FIGURE l8. Variation of Second-Order Exciting 
Force Amplitudes with Depth, 
Semi-Circular Hull in Heave 
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that ’ and ' are most sensitive to finite depth at high 
d s 

values of 6 or in very shallow water (h/d £ 1.5). 

(2 ) 

Figure 19 shows the corresponding effect on y of 

finite depth for the same range of h/d. Note that at 

( 2 ) 

h/d = 1.4 6 has negligible effect on y^ \ 

Also the effect of decreasing depth causes a monotonic 
( 2 ) 

variation of y at fixed 6. 




FIGURE 19. Variation of Second-Order Exciting Force Phase 
Angles with Depth, Semi-Circular Hull in Heave 
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Figure 20 presents the ratio of to 

plotted vs 6 for various values of h/d. The curves peak 
close to 6 = 1 for all values of h/d studied. The curves 
peak because the first-order force P^^^ passes through a 
minimum very near 6=1 for all values of h/d studied. 
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FIGURE 20. 



Ratio of Second-Order Exciting Force Amplitude 
to First-Order Exciting Force Amplitude, 
Semi-Circular Hull in Heave 



A quantitative illustration of the physical signif- 
icance of the nonlinear effects (second-order) is readily 
made using figure 20. Assume a ship with an approximately 
circular hull form is being excited such that the maximum 
acceleration of the motion is 0.2 g's. The maximum time- 
dependent amplitude of the force is given by 

l^^^^lmax " , (130) 



which may be written 



|F(t) 



e f 



*max _ T . _ 

^ " iOT 



(131) 



Assume h/d = 1.5 and 6 = 1.0, the maximum ship acceleration 
(in g's) may be written 

= e S . (132) 

g 



Equation 132 implies e = 0.2 under the assumed conditions. 
Using figure 20 and equations 131 and 132 



|F(t) 



e f 



max _ T 

HD ^ 



+ (0.2)(1.75) = 1.35 



(133) 



Therefore, neglecting the second-order contribution could 
result in a thirty-five percent error in |F(t)|^^^ . 

Figure 21 further demonstrates the nonlinear (second- 
order) effects on F(t). The function F(t)/e f^^^ is plotted 
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vs t. The ship's displacement as a function of time is also 
shown on the plot for reference (cosine function of unit 
amplitude). The parameters chosen for the plot are h/d = 1.5, 
6 =1.0 and e = 0.2 (corresponds to an equivalent acceleration 
of 0.2 g’s). 




time 



FIGURE 21. Total Exciting Force vs. Time 



Figure 22 Illustrates, for the same values of the 
parameters, the nonlinear effects on the gravity wave 
generated by the ship's disturbance. The figure presents 
the sum of the dimensionless first-order wave amplitude and 
the second-order Stokes' wave (both have celerity c^). 
Separately, the second-order wave from (celerity C 2 ) is 

plotted. Both curves are plotted vs x/b for fixed time. 




FIGURE 22. Second-Order V/ave Form Effects 
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VI. CONCLUSIONS 



It is concluded that this work has demonstrated the 
feasibility of using the finite element method to solve the 
linear and nonlinear problems of the type presented. Further, 
the finite element method is shown to be quite versatile in 
that various boundary conditions may be applied to the 
boundary value problems studied, with only minor changes in 
the computational scheme. In addition, the Isoparametric 
elements chosen are capable of representing arbitrary hull 
shapes . 

It is further concluded that a first-order study of the 
effect of finite depth on the ship motions studied demon- 
strates that little change in the solutions occur. However, 
it is asserted that the second-order results obtained 
demonstrate that nonlinear (second-order) effects on heaving, 
motions are most significant in very shallow water; particu- 
larly for motions at a frequency corresponding very closely 
to the natural frequency of free oscillations of the floating 
body. 

The present state of digital computer capacity appears to 
preclude the extension of the present problem to three- 
dimensions. However, the extension of this study to other 
modes of motion (sway and roll) is possible as demonstrated 
for the linear sway problem. Further, the possibility of 
studying non-symmetric , two-dimensional hull forms is 
suggested . 
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An adequate finite element mesh to represent the region 
for the nonlinear problems studied tends to tax the core 
capacity of the computer used. Nevertheless, meaningful 
results were obtained through careful design of the computa- 
tional scheme. 
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APPENDIX A 



Complex Algebra 



The development which follows determines the correct 
form of the complex numbers w^ and W 2 In the following 
expression 



Re{z^e^^^}Re{z 2 e^‘^^} = Re{w^ + W 2 e^^^^} 



(Al) 



z^ and Z 2 are complex numbers with moduli r^ and r 2 and 
arguments 9^ and 62 . It follows using appropriate trigono- 
metric identies that 

Re{r3,e“lel<’‘}Re{r2e“2ei°t} = 

l'*'COs 2 at ) + sin 9 ^sin 02 Js(l-cos 2 at) 

•p p 

- sin( 6 ^+ 02 )%sin 2 at] = [cos( 9 j^- 02 ) + cos ( 9 ^+ 02 + 2 at ) ] (A2) 



Equation A2 implies that 



w. 



Z1Z2 



W 2 = 



Z1Z2 



(A3) 



The symbol ('") implies conjugation. 

Equation A3 shows the origin of the factor % in various 
places in the text. 
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APPENDIX B 



Isoparametric Elements 






i 



FIGURE B1 . Element Mapping from C,ri Plane to x,y Plane 



The essential idea in the isoparametric element lies 
in the following equations which define a map from the 
plane to the x,y plane. 



N®x® , 


(Bl) 


N®y® 


(B2) 



The vector N® is a 1 x 12 row vector of element shape functions 
one for each node in figure Bl. The vectors x® and y® are 



\ 
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12 X 1 column vectors representing the respective x and y 
coordinates of the nodes In the x,y plane (figure Bl). 

The shape functions are conveniently defined using the 
conventions of Zienkiewicz [^5]. Let 

, Pq = r\T)^ , (B3) 

where C and n are the coordinates of any point of the square 
element In the Cjh plane. and are the coordinates of 

the 1^ node. The shape functions are given by the following 
set of equations. 

Corner Nodes: 

N® = ^ (1 + Co)(l + n^^)[-10 + + n^)] , (B^^) 

Edge Nodes: 

for q = ±1, , 

N® = ^ (1 + Co)(l - n^)(l + 9r)o) , (B5) 

for > hjL = ±1 > 

N® = ^ (1 + riQ)(l - 5^)(1 + 9 Cq) (B6) 
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The element Is defined to be isoparametric when the same 
shape functions are used to achieve the element mapping as 
are used to define the field function in the x,y plane. 

There are certain consistency requirements on the element 
shape functions as discussed by Zienklewicz [^5]. The shape 
functions chosen satisfy them all. 

The mapping has the advantage of being able to represent 
a region with curvilinear sides. The area and line integrals 
given in part III are calculated by the expressions given 
below. 

Prom the calculus, the following relation is known to hold 





r 




N® 






= J 


-X 


/ 


N® 




N® 








1 -y \ 



(B7) 



J is the real 2x2 Jacobian matrix of the transformation. 



riS e 
N^x 



N X 
~n~ 



ri© © 

N^y 



IvT© © 

N y 



(B8) 



It follows, provided J is not singular, that 







N® 


~x 


= j-i 




N® 




N® 


_~y _ 







(B9) 
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Then, for example, equation 86 in part III may be written at 
the element level as 



T T 

= / [N^ ] 

= Re ~ 





dRe 



Using equations B8, B9, and 



(BIO) 



dRe = dxdy = |J| dCdn , 



(BID 



equation B9 may be written 



+1 +1 

H® = / / [N® 

~ -1 -1 



N®*^] [J"^] 

z z 



N® 

~n 



|J| d^dn 



(B12) 



Equation B12 demonstrates the method of the evaluation 
of the matrix H® or In general the matrix H. The Integrations 
are accomplished using slx-polnt Gauss quadrature. 



APPENDIX C 



Computer Program Economization 




FIGURE Cl. Schematic of Coefficient Matrix 

The economizations in the computer solution of the linear 
system of equations defined by equation 102 in part IV are 
achieved by observing certain properties of the matrix 
First all interior nodes of a given mesh contribute only real 
numbers to This is also true of nodes on S^, S2 and 

Only the nodes on contribute complex entries in . 

Further only nodes on and contribute entries which are 
frequency dependent. 

The nodes on and typically comprise no more than 
fifteen percent of Further, the nodes on alone 
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comprise nor more than five percent of the total rows in . 

Therefore if is rearranged with equations for nodes on 

at the bottom, equations for nodes directly above the 
S|^ equations, and all other rows of the matrix filling the 
remaining part of , the resulting matrix would appear as 

schematically represented in figure Cl. The region designated 
with a Roman numeral I (uncrosshatched) indicates the portion 
of the matrix completely real and frequency independent. 

Region II indicates the portion of the matrix which has 

elements that are frequency dependent. And, Region III 

indicates the portion of the matrix which is complex and 

frequency dependent. If the matrix is formed as shown 

in figure Cl, eighty-five percent (approximately) of it may 

be eliminated only once using real arithmetic. This is 

2 

possible because the matrix -a Qq does not have to be added 
to- Region II until after Region I has been reduced by Gauss 
elimination, and the same fact applies to the matrix 
(ia/c^)D in Region III. Then, after Region I is eliminated, 
the matrix -o Qq is added to Region II (after saving the 
previous results) and is reduced to upper triangular 

form through Region II using real arithmetic. Now, if this 
result is saved and then the matrix (ia/c^)D is added to 
Region III only a small set of complex equations (usually 
about 25) is left to solve. 

After the final elimination, the resulting matrix is in 
complete upper-triangular form and ready for a back- 
substitution process. The right-hand side vector is complex 
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and frequency dependent In general. However, the required 
Information to perform the forward substitution process on 
it is contained in the final form of (upper triangular) . 

Therefore, the forward substitution on Ibar^^^ may be done 
after is processed. Then a back-substitution of the 

resulting system using mixed-mode computer arithmetic yields 
the final solution vector . 

A close examination of the process just described reveals 
that for a whole family of frequencies only a small portion 
of has to be reassembled and eliminated again. These 

matrix properties contribute toward a sizeable saving in 
computer time. 
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